function [owner_profit1, owner_profit2, owner_profit, ...
    mean_p1, mean_p2, mean_p, ttq1, ttq2, ttq, ...
    cons1, cons2, cons, nprod1, nprod, is_same, ...
    prod_sim_fc, added_prod, dropped_prod, ind_cf]...
    = prod_equi_sim(config_init_by_owner, ...
    prod_id_by_owner, prod_attributes_m, brewer_name_by_owner, merge_owner, ...
    fixed_prod_ind_id, fixed_prod_entry_m, change_prod, prod_sim_fc_all, prod_sim_fc0, n_fc_sim, pop, ...
    para_fc0, para_profit, para_in_state, para_first, para_craft, para_mkt_sig, cost_reduction, merged_craft, ...
    list_owner, sim_pcoeff, setup, ...
    para_nonlinear, randomdraw, ...
    rand_income, income_norm, starting_point_option, ...
    owner_seq, check_sim_fc, mktid, np)


%% setup

setup.BR_method = 'one_product'; 
setup.unique_ind_ijk = [];
setup.recover_ind_ijk = [];
setup.apple = [];

n_owner = length(prod_id_by_owner);
max_iter_count = 5;

para_first_normalized = para_first/para_profit/pop;% in the unit of per-capita $

for no = 1 : n_owner
    if isempty(config_init_by_owner{no})
        config_init_by_owner{no} = false(0,1);
    end
end

%% simulate fixed cost shocks consistent with the data being an eqm
tic; fprintf(1, 'np=%1.0f, mktid=%1.0f, start with FC draws\n', np, mktid);
display_setup.display = false; 

owner_var_profit = prod_config_profit(prod_attributes_m, ...
    fixed_prod_ind_id, fixed_prod_entry_m, config_init_by_owner, prod_id_by_owner, ...
    sim_pcoeff, setup, ...
    para_nonlinear, randomdraw, ...
    rand_income, income_norm);

if all(owner_var_profit == 0) || any(isnan(owner_var_profit)); error('this market should have endogenous products'); end

%(1) range of fc in per-capita $
fc_range = cell(n_owner, 1);
fixed2 = cell(n_owner, 1);
fixed2_entry = cell(n_owner, 1);

for no = 1 : n_owner
    n_prod_no = length(config_init_by_owner{no});
    
    fc_range_no = nan(n_prod_no, 2);
    for j = 1 : n_prod_no
        config_by_owner_dev = config_init_by_owner;
        config_by_owner_dev{no}(j) = 1 - config_init_by_owner{no}(j);

        owner_var_profit_dev = prod_config_profit(prod_attributes_m, ...
            fixed_prod_ind_id, fixed_prod_entry_m, config_by_owner_dev, prod_id_by_owner, ...
            sim_pcoeff, setup, ...
            para_nonlinear, randomdraw, ...
            rand_income, income_norm);

        profit_diff = owner_var_profit_dev(no) - owner_var_profit(no);

        if config_init_by_owner{no}(j)==0 % lower bound of FC
            fc_range_no(j, :) = [profit_diff, profit_diff+1]; % add $1/mkt/product to ensure bounds are valid
        elseif config_init_by_owner{no}(j)==1 % upper bound of FC
            fc_range_no(j, :) = [-profit_diff-1, -profit_diff]; % subtract $1/mkt/product to ensure bounds are valid
        else
            error('error: binary indexing');
        end
    end

    fc_range{no} = fc_range_no;
end

if ~isempty(prod_sim_fc0) % if a market with no products is included
    prod_sim_fc = prod_sim_fc0;
else
    prod_sim_fc = cell(n_owner, 1);

    for no = 1 : n_owner
        % (2) range of fc unobservable in the unit of "fc shock std deviation"
        %  get the covariates
        n_prod_no = length(config_init_by_owner{no});
        
        in_state_no = nan(n_prod_no, 1);
        
        craft_no = nan(n_prod_no, 1);

        list_brewer_no = unique(brewer_name_by_owner{no});
        n_brewer_no = length(list_brewer_no);

        if n_brewer_no ==1
            tmp = sum(config_init_by_owner{no}) + 1 - config_init_by_owner{no};
            isfirst_no = (tmp==1); %for a product in the portfolio, isfirst_no = 1 if dropping it leads to droping the scope parameter (i.e., sum(config_init_by_owner{no})=1);
                    %for a product not in the portfolio, isfirst_no = 1 if adding it leads to adding the scope parameter (i.e.,sum(config_init_by_owner{no})=0);
        else
            tmp = nan(n_prod_no, 1);
            for nb =1 : n_brewer_no
                ind_nb = ismember(brewer_name_by_owner{no}, list_brewer_no(nb));

                tmp(ind_nb) = sum(config_init_by_owner{no}(ind_nb)) + ...
                    1 - config_init_by_owner{no}(ind_nb);
            end
            isfirst_no = (tmp==1);
        end

        for j = 1 : n_prod_no
            in_state_no(j) = prod_attributes_m{prod_id_by_owner{no}(j)}.in_state;

            if isequal(prod_attributes_m{prod_id_by_owner{no}(j)}.brewer, 'goose island beer co')
                craft_no(j) = 0;
            else
                craft_no(j) = prod_attributes_m{prod_id_by_owner{no}(j)}.craft;
            end
        end

        % the mean portion of the fixed cost in the unit of "fc shock std deviation", excluding the scope parameter
        % note: the fixed cost function in the code will account the scope parameter 
        if para_mkt_sig == 0
            mean_fc_no = para_fc0 + in_state_no*para_in_state + craft_no*para_craft;
        else
            mean_fc_no = para_fc0 + in_state_no*para_in_state + craft_no*para_craft + para_mkt_sig * randn(1,1);
        end

        % the range of fc unobservable
        fc_unobs_range_no = fc_range{no}*pop*para_profit - mean_fc_no - isfirst_no*para_first;

        % fix the entry outcomes of some products
        fixed2{no} = (fc_unobs_range_no(:, 1).*fc_unobs_range_no(:, 2)>0)... 
            &(min(abs(fc_unobs_range_no), [], 2)>3);

        fixed2_entry{no} = config_init_by_owner{no}(fixed2{no});

        %(3) simulate fc draws (in the unit of per-capita $)
        prod_sim_fc_no = nan(n_prod_no, n_fc_sim);        

        if check_sim_fc % check if there is profitable deviation
            k = 1;
            ns = 1;
            firm_seq_for_checking = no;
            while k<=size(prod_sim_fc_all{no}, 2) && ns<=n_fc_sim
                sim_ns = simtnorm_double(fc_unobs_range_no(:, 1), fc_unobs_range_no(:, 2), prod_sim_fc_all{no}(:, k));
                fc_sim_ns = (sim_ns + mean_fc_no)/para_profit/pop;% in the unit of per-capita $

                prod_sim_fc_for_checking = cell(n_owner, 1);
                prod_sim_fc_for_checking{no} = fc_sim_ns;

                nfc_for_checking = 1;

                %compute the best response based on a starting point different from the configuration at the data
                % and check whether the computed best response is the same as the configuration at the data
                config_init_for_checking = config_init_by_owner;
                config_init_for_checking{no} = true(size(config_init_by_owner{no}));
        
                config_binary_nb = prod_equi_swap(...
                    config_init_for_checking, prod_attributes_m,...
                    fixed_prod_ind_id, fixed_prod_entry_m, fixed2, fixed2_entry, ...
                    prod_sim_fc_for_checking, para_first_normalized, nfc_for_checking, ...
                    firm_seq_for_checking, sim_pcoeff,   ...
                    max_iter_count, setup, ...
                    prod_id_by_owner, brewer_name_by_owner, ...
                    para_nonlinear, randomdraw, ...
                    rand_income, income_norm, [], display_setup);

                if ~isequal(config_binary_nb, config_init_by_owner) % adjust fc smln draws

                    adjust_shock_iter = 0;
                    while ~isequal(config_binary_nb, config_init_by_owner)
                        adjust_shock_iter = adjust_shock_iter + 1;

                        % adjust product by product
                        for kk = 1 : size(fc_sim_ns)
                            if adjust_shock_iter<1
                                if config_binary_nb{no}(kk)==1 && config_init_by_owner{no}(kk)==0
                                    sim_ns(kk) = simtnorm_double(sim_ns(kk), sim_ns(kk)+10, prod_sim_fc_all{no}(kk, k));
                                elseif config_binary_nb{no}(kk)==0 && config_init_by_owner{no}(kk)==1
                                    sim_ns(kk) = simtnorm_double(sim_ns(kk)-10, sim_ns(kk), prod_sim_fc_all{no}(kk, k));
                                end
                            else
                                % likely outlier; push to the boundary (effectively fix the product)
                                if config_binary_nb{no}(kk)==1 && config_init_by_owner{no}(kk)==0
                                    sim_ns(kk) = sim_ns(kk) + 100000;
                                elseif config_binary_nb{no}(kk)==0 && config_init_by_owner{no}(kk)==1
                                    sim_ns(kk) = sim_ns(kk) - 100000;
                                end
                            end
                        end
                        fc_sim_ns = (sim_ns + mean_fc_no)/para_profit/pop;% in the unit of per-capita $
                        prod_sim_fc_for_checking{no} = fc_sim_ns;
                        config_binary_nb = prod_equi_swap(...
                            ...
                            config_init_for_checking, prod_attributes_m,...
                            fixed_prod_ind_id, fixed_prod_entry_m, fixed2, fixed2_entry, ...
                            prod_sim_fc_for_checking, para_first_normalized, nfc_for_checking, ...
                            firm_seq_for_checking, sim_pcoeff,   ...
                            max_iter_count, setup, ...
                            prod_id_by_owner, brewer_name_by_owner, ...
                            para_nonlinear, randomdraw, ...
                            rand_income, income_norm, [], display_setup);
                    end

                    if ismember(no, merged_craft)
                        prod_sim_fc_no(:, ns)  = max(0, fc_sim_ns - cost_reduction/para_profit/pop);% in the unit of per-capita $
                    else
                        prod_sim_fc_no(:, ns) = fc_sim_ns;
                    end

                    ns = ns + 1;
                    k = k + 1;                    
                else
                    if ismember(no, merged_craft)
                        prod_sim_fc_no(:, ns)  = max(0, fc_sim_ns - cost_reduction/para_profit/pop);% in the unit of per-capita $
                    else
                        prod_sim_fc_no(:, ns) = fc_sim_ns;
                    end

                    ns = ns + 1;
                    k = k + 1;
                end

            end            
        else

            if ismember(no, merged_craft)
                for ns = 1 : n_fc_sim
                    sim_ns = simtnorm_double(fc_unobs_range_no(:, 1), fc_unobs_range_no(:, 2), prod_sim_fc_all{no}(:, ns));
                    fc_sim_ns = max(0, (sim_ns + mean_fc_no - cost_reduction)/para_profit)/pop;% in the unit of per-capita $
                    prod_sim_fc_no(:, ns) = fc_sim_ns;
                end
            else
                for ns = 1 : n_fc_sim
                    sim_ns = simtnorm_double(fc_unobs_range_no(:, 1), fc_unobs_range_no(:, 2), prod_sim_fc_all{no}(:, ns));
                    fc_sim_ns = (sim_ns + mean_fc_no)/para_profit/pop;% in the unit of per-capita $
                    prod_sim_fc_no(:, ns) = fc_sim_ns;
                end
            end

        end
        prod_sim_fc{no} = prod_sim_fc_no;
    end
end
fprintf(1, 'np=%1.0f, mktid=%1.0f, done with FC draws, time = %1.2f\n', np, mktid, toc);

%% counterfactual setup: merger
tic; fprintf(1, 'np=%1.0f, mktid=%1.0f, start with CF computation\n', np, mktid);

prod_id_by_owner_cf = prod_id_by_owner;
config_init_by_owner_cf = config_init_by_owner;
fixed2_cf = fixed2;
fixed2_entry_cf = fixed2_entry;
brewer_name_by_owner_cf = brewer_name_by_owner;
prod_attributes_cf = prod_attributes_m;
prod_sim_fc_cf = prod_sim_fc;
owner_seq_cf = owner_seq;

if length(merge_owner)>1 
    % use the first firm to represent the merging firms
    merge_owner_id = find(ismember(list_owner, merge_owner));
    
    for nm = 2 : length(merge_owner)
        prod_id_by_owner_cf{merge_owner_id(1)} = ...
            [prod_id_by_owner_cf{merge_owner_id(1)}; prod_id_by_owner{merge_owner_id(nm)}];
        prod_id_by_owner_cf{merge_owner_id(nm)} = [];

        config_init_by_owner_cf{merge_owner_id(1)} = ...
            [config_init_by_owner_cf{merge_owner_id(1)}; config_init_by_owner{merge_owner_id(nm)}]>0;
        config_init_by_owner_cf{merge_owner_id(nm)} = true(0,1);

        fixed2_cf{merge_owner_id(1)} = ...
            [fixed2_cf{merge_owner_id(1)}; fixed2{merge_owner_id(nm)}]>0;
        fixed2_cf{merge_owner_id(nm)} = true(0, 1);

        fixed2_entry_cf{merge_owner_id(1)} = ...
            [fixed2_entry_cf{merge_owner_id(1)}; fixed2_entry{merge_owner_id(nm)}]>0;
        fixed2_entry_cf{merge_owner_id(nm)} = true(0, 1);

        brewer_name_by_owner_cf{merge_owner_id(1)} = ...
            [brewer_name_by_owner_cf{merge_owner_id(1)}; brewer_name_by_owner{merge_owner_id(nm)}];
        brewer_name_by_owner_cf{merge_owner_id(nm)} = [];

        prod_sim_fc_cf{merge_owner_id(1), :} = ...
            [prod_sim_fc_cf{merge_owner_id(1)}; prod_sim_fc{merge_owner_id(nm), :}];
        prod_sim_fc_cf{merge_owner_id(nm)} = [];
    end

    % change owner_ind in prod attributes
    for j = 1 : length(prod_attributes_m)
        if sum(prod_attributes_m{j}.owner_ind(merge_owner_id(2:end)))>0
            prod_attributes_cf{j}.owner_ind(merge_owner_id(2:end)) = 0;
            prod_attributes_cf{j}.owner_ind(merge_owner_id(1)) = 1;
        end
    end
else
    merge_owner_id = [];
end


%% merger/divestiture simulation
%(1) at the data
[~, ~, ~, owner_ind_n1, ~, ~, ~, share1,...
    cons_surplus1, ~, ~, carrier_profit1, p_sim1, ...
    w_sim1, share_sim1, prod_var_profit, brand_in_mkt, ~, craft1] = prod_config_profit(prod_attributes_m, ...
    fixed_prod_ind_id, fixed_prod_entry_m, config_init_by_owner, prod_id_by_owner, ...
    sim_pcoeff, setup, ...
    para_nonlinear, randomdraw, ...
    rand_income, income_norm);

owner_profit1 = nan(n_fc_sim, 2); %col1: producer surplus, col2: merging firms' profit

if isempty(merge_owner_id)
    ind_merged_craft_brands = false(size(owner_ind_n1, 1), 1);
else
    ind_merged_craft_brands = sum(owner_ind_n1(:, merge_owner_id), 2)>0;
end

for n = 1 : n_fc_sim
    [~, fc_prod_n1] = ...
        compute_fc(length(prod_attributes_m), config_init_by_owner, ...
        brewer_name_by_owner, prod_id_by_owner, fixed_prod_ind_id, ...
        fixed_prod_entry_m, para_first_normalized, prod_sim_fc, n);
    fc_prod_n1 = max(fc_prod_n1, 0);
    prod_profit = prod_var_profit - fc_prod_n1;

    owner_profit1(n, :) = pop*[sum(prod_profit), sum(prod_profit(ind_merged_craft_brands))];
end

nprod1 = [size(craft1, 1), ...
    sum(craft1), ...
    sum(craft1(ind_merged_craft_brands))];

mean_p1 = [mean(sum(p_sim1.*share_sim1, 1)./sum(share_sim1, 1)), ...
    mean(sum(p_sim1(craft1>0, :).*share_sim1(craft1>0, :), 1)....
    ./sum(share_sim1(craft1>0, :), 1)), ...
    mean(sum(p_sim1(ind_merged_craft_brands & craft1>0, :).*share_sim1(ind_merged_craft_brands & craft1>0, :), 1)....
    ./sum(share_sim1(ind_merged_craft_brands & craft1>0, :), 1))];

ttq1 = round(pop*[sum(share1), ...
    sum(share1(craft1>0)), ...
    sum(share1(ind_merged_craft_brands & craft1>0))]);

cons1 = cons_surplus1*pop;

[~, ~, ~, owner_ind_n2, ~, ~, ~, share2,...
    cons_surplus2, ~, ~, carrier_profit2, p_sim2, ...
    w_sim2, share_sim2, prod_var_profit2, brand_in_mkt2, ~, craft2] = prod_config_profit(prod_attributes_cf, ...
    fixed_prod_ind_id, fixed_prod_entry_m, config_init_by_owner_cf, prod_id_by_owner_cf, ...
    sim_pcoeff, setup, ...
    para_nonlinear,randomdraw,  ...
    rand_income, income_norm);

owner_profit2 = nan(n_fc_sim, 2);
for n = 1 : n_fc_sim
    [~, fc_prod_n2] = ...
        compute_fc(length(prod_attributes_cf), config_init_by_owner_cf, ...
        brewer_name_by_owner_cf, prod_id_by_owner_cf, fixed_prod_ind_id, ...
        fixed_prod_entry_m, para_first_normalized, prod_sim_fc_cf, n);
    fc_prod_n2 = max(fc_prod_n2, 0);
    prod_profit2 = prod_var_profit2 - fc_prod_n2;
    owner_profit2(n, :) = pop*[sum(prod_profit2), sum(prod_profit2(ind_merged_craft_brands))];
end

mean_p2 = [mean(sum(p_sim2.*share_sim2, 1)./sum(share_sim2, 1)), ...
    mean(sum(p_sim2(craft2>0, :).*share_sim2(craft2>0, :), 1)....
    ./sum(share_sim2(craft2>0, :), 1)), ...
    mean(sum(p_sim2(ind_merged_craft_brands & craft2>0, :).*share_sim2(ind_merged_craft_brands & craft2>0, :), 1)....
    ./sum(share_sim2(ind_merged_craft_brands & craft2>0, :), 1))];

ttq2 = round(pop*[sum(share2), ...
    sum(share2(craft2>0)), ...
    sum(share2(ind_merged_craft_brands & craft2>0))]);

cons2 = cons_surplus2*pop;

%(3) full counterfactual with product changes
ind_cf = cell(n_fc_sim, 1);
if change_prod
    for nfc_cf = 1 : n_fc_sim
        display_setup.display = true; 
        display_setup.np = n;
        display_setup.mktid = mktid;
        display_setup.nfc_cf = nfc_cf;

        ind_cf{nfc_cf} = prod_equi_swap(...
            ...
            config_init_by_owner_cf, prod_attributes_cf,...
            fixed_prod_ind_id, fixed_prod_entry_m, fixed2_cf, fixed2_entry_cf, ...
            prod_sim_fc_cf, para_first_normalized, nfc_cf, ...
            owner_seq_cf, sim_pcoeff,   ...
            max_iter_count, setup, ...
            prod_id_by_owner_cf, brewer_name_by_owner_cf, ...
            para_nonlinear, randomdraw, ...
            rand_income, income_norm, starting_point_option, display_setup);
    end
else
    for nfc_cf = 1 : n_fc_sim
        ind_cf{nfc_cf} = config_init_by_owner_cf;
    end
end

is_same = nan(n_fc_sim, 1);
for nfc_cf = 1 : n_fc_sim
    is_same(nfc_cf) = isequal(ind_cf{nfc_cf}, config_init_by_owner_cf);
end

mean_p = nan(n_fc_sim, 3);
ttq = nan(n_fc_sim, 3);
cons = nan(n_fc_sim, 1);
owner_profit = nan(n_fc_sim, 2);
nprod = nan(n_fc_sim, 3);
added_prod = cell(n_fc_sim, 1);
dropped_prod = cell(n_fc_sim, 1);
owner_ind = cell(n_fc_sim, 1);

for n = 1 : n_fc_sim
    config_dn = ind_cf{n};
    [~, ~, ~, owner_ind_n, ~, ~, ~, share,...
        cons_surplus, ~, ~, ~, p_sim, ...
        ~, share_sim, prod_var_profit_n, brand_in_mkt_nfc, ~, craft_nfc] = prod_config_profit(prod_attributes_cf, ...
        fixed_prod_ind_id, fixed_prod_entry_m, config_dn, prod_id_by_owner_cf, ...
        sim_pcoeff, setup, ...
        para_nonlinear, randomdraw, ...
        rand_income, income_norm);

    owner_ind{n} = owner_ind_n;

    ind_merged_craft_brands_nfc = sum(owner_ind_n(:, merge_owner_id), 2)>0;
    
    added_prod{n} = setdiff(brand_in_mkt_nfc, brand_in_mkt);
    
    dropped_prod{n} = setdiff(brand_in_mkt, brand_in_mkt_nfc);

    nprod(n, :) = [size(craft_nfc, 1), sum(craft_nfc),...
        sum(ind_merged_craft_brands_nfc & craft_nfc>0)];

    mean_p(n, :) = [mean(sum(p_sim.*share_sim, 1)./sum(share_sim, 1)), ...
        mean(sum(p_sim(craft_nfc>0, :).*share_sim(craft_nfc>0, :), 1)....
        ./sum(share_sim(craft_nfc>0, :), 1)), ...
        mean(sum(p_sim(ind_merged_craft_brands_nfc & craft_nfc>0, :).*share_sim(ind_merged_craft_brands_nfc & craft_nfc>0, :), 1)....
        ./sum(share_sim(ind_merged_craft_brands_nfc & craft_nfc>0, :), 1))];

    ttq(n, :) = round(pop*[sum(share), ...
        sum(share(craft_nfc>0)), ...
        sum(share(ind_merged_craft_brands_nfc & craft_nfc>0))]);

    cons(n) = cons_surplus*pop;

    [~, fc_prod] = ...
        compute_fc(length(prod_attributes_m), config_dn, ...
        brewer_name_by_owner_cf, prod_id_by_owner_cf, fixed_prod_ind_id, ...
        fixed_prod_entry_m, para_first_normalized, prod_sim_fc_cf, n);

    fc_prod = max(0, fc_prod);
    prod_profit_n = prod_var_profit_n - fc_prod;
    owner_profit(n, :) = pop*[sum(prod_profit_n), sum(prod_profit_n(ind_merged_craft_brands_nfc))];
end


%(4) counterfactual with no entry
firm_id_all = find(sum(owner_ind_n1, 1));
merge_owner_id_in_mkt = intersect(merge_owner_id, firm_id_all);
non_merge_firm_id = setdiff(firm_id_all, merge_owner_id);

craft_firm_id_all = find(sum(owner_ind_n1(craft1, :), 1));
non_craft_firm_id_all = find(sum(owner_ind_n1(craft1, :), 1)==0);
craft_firm_id_merge = intersect(craft_firm_id_all, merge_owner_id_in_mkt);

n_all = length(firm_id_all);
n_craft = length(craft_firm_id_all);
n_firm = [n_all, n_craft];

non_merge_firm_id = setdiff(firm_id_all, merge_owner_id);
n_non_merge_firm_id = length(non_merge_firm_id);

n_firm_sim = nan(n_fc_sim, 1);
added_prod_by_non_merge = zeros(n_fc_sim, 1);
added_prod_by_non_merge_entry = zeros(n_fc_sim, 1);
ind_cf_no_de_novo = cell(n_fc_sim, 1);

for n = 1 : n_fc_sim
    ind_cf_no_de_novo{n} = ind_cf{n};
    owner_ind_n = owner_ind{n};
    
    non_merge_firm_id_cf = setdiff(find(sum(owner_ind_n, 1)), merge_owner_id);
    n_non_merge_firm_id_n = length(non_merge_firm_id_cf);

    n_firm_sim(n, :) = n_non_merge_firm_id_n-n_non_merge_firm_id;
    for k = non_merge_firm_id_cf
        added_prod_by_non_merge(n) = added_prod_by_non_merge(n) + ...
            double((sum(config_init_by_owner{k})>0)*(sum(ind_cf{n}{k}) - sum(config_init_by_owner{k})));
        added_prod_by_non_merge_entry(n) = added_prod_by_non_merge_entry(n) + ...
            double((sum(config_init_by_owner{k})==0)*(sum(ind_cf{n}{k}) - sum(config_init_by_owner{k})));
        if (sum(config_init_by_owner{k})==0)*(sum(ind_cf{n}{k}) - sum(config_init_by_owner{k}))>0
            ind_cf_no_de_novo{n}{k} = false(size(ind_cf_no_de_novo{n}{k}));
        end
    end
end

prod_owner_ind_no_craft = config_init_by_owner;

for n_att = 1 : length(config_init_by_owner)
    prod_owner_ind_no_craft{n_att} = zeros(size(config_init_by_owner));
end

% save surplus
[~, ~, ~, ~, ~, ~, ~, ~,...
    cons_surplus_no_craft] = prod_config_profit(prod_attributes_m, ...
    fixed_prod_ind_id, fixed_prod_entry_m, prod_owner_ind_no_craft, prod_id_by_owner, ...
    sim_pcoeff, setup, ...
    para_nonlinear, randomdraw, ...
    rand_income, income_norm);

cons_no_de_novo = nan(n_fc_sim, 1);
for n = 1 : n_fc_sim
    [~, ~, ~, ~, ~, ~, ~, ~,...
        cons_n_nd] = prod_config_profit(prod_attributes_cf, ...
        fixed_prod_ind_id, fixed_prod_entry_m, ind_cf_no_de_novo{n}, prod_id_by_owner_cf, ...
        sim_pcoeff, setup, ...
        para_nonlinear, randomdraw, ...
        rand_income, income_norm);

    cons_no_de_novo(n) = cons_n_nd;
end

fprintf(1, 'np=%1.0f, mktid=%1.0f, done with CF computation, time = %1.2f\n', np, mktid, toc);

list_variables = {'n_owner','ind_cf','config_init_by_owner','prod_id_by_owner', 'prod_attributes_m', 'para_fc0', 'para_in_state', 'para_craft', 'para_profit', 'merge_owner_id', 'owner_profit1', 'owner_profit2', 'owner_profit', 'mean_p1', 'mean_p2', 'mean_p', 'ttq1', 'ttq2', 'ttq', 'cons1', 'cons2', 'cons', 'nprod1', 'nprod', 'n_firm', 'n_firm_sim', 'merge_owner_id_in_mkt', 'added_prod', 'dropped_prod', 'ind_cf', 'pop', 'cons_surplus_no_craft', 'cons_no_de_novo', 'added_prod_by_non_merge_entry', 'added_prod_by_non_merge'};
if isempty(merge_owner)
    save(['cf_results_', num2str(mktid), '_', num2str(np), '_pre'], list_variables{:});
else
    save(['cf_results_', num2str(mktid), '_', num2str(np), '_post'], list_variables{:});
end
